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ABSTRACT 


Results are given of a comparative study of numerical procedures for 
computing solid wall boundary points in supersonic inviscid flow calculations. 
Twenty five different calculation procedures were tested on two sample problems, 
a simple expansion wave and a simple compression (two dimensional steady flow) . 
A new simple but accurate calculation procedure was developed. The merits and 
shortcomings of the various procedures are thoroughly discussed along with 
complications for three dimensional and time dependant flows. 
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SECTION 1 


INTRODUCTION 

The purpose of this study is twofold: (1) to review and evaluate tech- 

niques for handling boundary points in the numerical solution of supersonic 
steady flow fields with finite difference procedures, and (2) to improve an 
existing or develop a new technique which is simple to implement and accurate. 

In order to clarify and simplify the study as much as possible, attention is 
primarily focused on two-dimensional, steady, supersonic flow of an ideal gas. 
Although some implicit procedures for handling boundary conditions are considered, 
the object is to consider methods that fit in well with explicit differencing 
procedures (particularly that of MacCormack) for the interior. 

The study is organized into three parts: (1) review, (2) evaluation, 

and (3) development. A preliminary review of current procedures was first per- 
formed, and it resulted in the categorization of the various procedures according 
to the basic principles incorporated. Then representative techniques from each 
category were compared on problems for which the exact solution is known. 

Unless otherwise specified, a two-dimensional steady flow is implied, with 
the flow direction globally from left to right. Axes and velocity components are 
(x,y) and (u,v) respectively, with the x axis positive to the right. Thus, a 
numerical solution is obtained by marching from the left to the right. 

This report begins (Section 2) with a brief review and categorization of 
computational procedures which are now in use, including one developed in this 
study. Then (Section 3) we discuss the most important elements of a computational 
procedure, with emphasis on boundary points. Comparative results, and a discussion 
of the significance, are given in Section 4. The report is concluded (Section 5) 
with an overview of the results of the study. 

It became evident early that the procedures currently in use which give 
acceptable accuracy tend, in general, to be cumbersome to implement in a computer 
code. Also, they often require relatively long computing times per mesh point 
when compared to interior points. The latter is not a strong objection, since it 
is in the boundary calculations that the particular solution is selected from the 
infinity of possible solutions to the partial differential equations. However, it 
is important in the sense that when boundary point calculations comprise a signifi- 
cant portion of the computing time, improving them can result in significant sav- 
ings in computational expense. 



Because of the obvious need for simpler, faster, accurate procedures, 
efforts were directed toward understanding the factors determining a good boundary 
point computation procedure, particularly in the context of a predictor /corrector 
scheme at interior points. Subsequently, a new scheme for computing boundary 
points was developed. This scheme is very simple, requires little computing time, 
and seems to be quite accurate. It is discussed in Sections 2, 3, and 4. 
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SECTION 2 
REVIEW 


2.1 CATEGORIZATION OF METHODS 

The various procedures have been classified in five categories: reflec- 
tion, explicit differencing (one-sided derivatives) , implicit, characteristics, 
and miscellaneous. Although the last category is somewhat of a grab-bag, it 
should not be viewed disparagingly, for that is probably where the most satis- 
factory procedure (from the standpoint of simplicity and accuracy) will be 
found. 


2.2 DISCUSSION OF VARIOUS PROCEDURES TO DATE 
2.2.1 Reflection 

The use of this procedure is basically physically motivated. The wall 
boundary of an inviscid flow is essentially a slip stream, and in this sense 
there is no basic difference between a streamline ?long a solid wall and the 
same streamline which is located within a fluid on one side of a slip stream. 
On this basis, one reasons that, for computational convenience, he can 
replace his original problem with its solid boundary with another problem 
having an extended flow field which (hopefully) will include a slip stream 
where there was originally a solid wall (sketch) . 




Original Flow Field 



Fictitiously Extended Flow Field 


Sketch 1 The Concept of Reflection Procedures 
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Knowing the solution, at a given initial line, somewhat into the extended flow,* ** 

t 

it is possible to continue the solution downstream accurately by treating the 
slip stream as a regular part of the interior flow. So the problem of computing 
the solution at the solid boundary has disappeared. Unfortunately, it h..s been 
replaced by the problem of determining the fictitiously extended flow field be- 
low the solid wall. 

In looking at two-dimensional problems with straight walas, it is not 
difficult to construct an extended flow field which will effectively simulate 
the solid boundary, ^nd there the idea of reflection enters naturally. In this 
case, at the wall the normal derivative of pressure vanishes along with the 
normal component of velocity. Therefore, the extended flow field can be con- 
structed by simply reflecting the actual flow field about the solid wall. All 
quantities are reflected evenly, except the normal velocity component, which 
suffers odd reflection in order to impose the requirement that it vanish at the 
solid wall. 

To carry this idea further, to flows with curved walls or to three- 

dimensional flows, leads to a real quagmire. Though there can be no argument 

in any particular case about the concept of an extended flow field, except in 

special cases the difficulty of constructing one which will yield an accurate 

solution to the posed problem is probably at least as difficult as it is to ob- 

** 

tain more directly a solution of comparable accuracy. 

It is evident from this discussion that the biggest shortcoming of re - 
flection procedures is that their accuracy depends so strongly on the local 
characteristics of the f 1 jw (and on the coordinate system used in the computation) , 
and no one has utilized procedures reflecting this fact with enough precis 
However, it should be recognized that the use of body oriented coordinate ;ms 

is an attempt to accomplish this. 

Reflection procedures have been used by a number of authors with varying 
degrees of success. The incorrect behavior of the sonic-line flow about an 
Apollo- like body as computed by Bohachevsky and Mates and Bohachevsky and 
Kostoff certainly results from poor solutions at boundary points.*** Reflection 
procedures in one form or another have been utilized by many other authors, in- 
cluding Bohachevasky and Rubin, Bur stein, Eaton, Kutler, Kutler and Lomax, 

* 

Note that this extended solution is not uniquely determined until one re- 
quires that the solution be continuous across the slipstream. That assump- 
tion is usually implicit. 

** 

This will become more evident in Section 4. 

*** 

Compare these solutions with those of Barnwell (1971) , for instance. 
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MacCormack, Serra, Tyler and Zumwalt, and Walker and Zumwalt. In the Latter 
report (pg 49) 'the following justification is made for utilizing a reflection 
procedure ... "The obvious advantage of the image point method (reflection) for 
treating boundary points is that it completely deletes the necessity of r’icta- 
ting a priori the v lues of derivatives at a surface." This statement is errone- 
ous and misleading nee exactly the opposite is true. It follows a statement 
(p. 42) that there .__*e "three ways by which one can take into account the presence 
of boundary when a finite difference method is employed." The latter statement 
is also erroneous (and was in 1966) since there are far more than three possibil- 
ities. Actually, it should be clear from this discussion that image point or 
reflection procedures themselves do impose normal derivatives at the boundary 
(usually = 0 for most quantities, ji 0 for normal velocity component). 

2.2.2 Explicit One-Sided Derivatives 

It has been implied above that it would be very convenient to solve 
boundary points exactly as interior points , and this is one of the driving 
forces behind the common utilization of reflection procedures. The problem is 
that virtually all satisfactory difference schemes for interior points require 
some form of centered differencing.* Howeve. , there is no reason to require 
that the idea of central differences be carried over to boundary points. In 
fact, it violates what we know to be the physical/mathematical structure cf 
the problem. At a solid wall, che solution along che wa2 ’ is determined solely 
by the interaction of the wall geometry and the interior flow. Therefore, though 
it is not quite as convenient, it is much more natural to think in terms of one- 
sided derivatives at the wall if one is embarked on a course of solving the 
partial differential equations along the wall . 

With this in mind, cne can easily construct procedures for integrating 
one or more of the differential equations along the wall with derivatives nor- 
mal to the wall in the computational space being computed hv one-sided differences. 
The actual procedures employed can vary somewhat, particularly with respect to 
how many and which equations are to be integrated in this manner. 

Such procedures have been employed by Barnwell (1970) , Grossman and 
Moretti, Li, Moretti (PIBAL 70-20), and Skoglund and Guy,** and they are discussed 
by Ciment, Eaton, and Walker and Zumwalt. Probably the most appealing aspect 


It need not be symmetric, and multi-step procedures are included in this 
rather broad interpretation of centered differencing. 

Though their procedur. is more than just utilizing one-sided derivatives. 
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of procedures based on simple one-sided differencing of the partial differential 

f 

equations is that, with the exception of reflection procedures, they most con- 
veniently fit into schemes for integrating the equations at interior points. 

2.2.3 Standard Predictor-Corrector Schemes - One-sided Derivatives 

A logical extention to the simplest explicit scheme with one-sided deriva- 
tives is to consider a predictor-corrector scheme with derivatives in both the 
predictor and corrector steps computed as one-sided differences on the same side. 
Such an approach is a logical extension of a predictor-corrector scheme at interior 
points tc boundary points. Though we can expect to have potentially more accurate 
calculations than those outlined in the preceding section, accuracy will be lim- 
ited by the fac that no direct coupling is made between the partial differential 
equations and the wall geometry. 

2.2.4 Implicit Schemes - One-sided Derivatives 

Implicit schemes, and mixed implicit/explicit schemes, are the same as 
corresponding completely explicit schemes except that the derivatives "normal" 
to the wall are computed with the values at the next x station (i.e., x = x Q + 

Ax). Thus, the solution depends on the "normal" derivatives which, in turn, 
depend on the solutior , and the solution must be determined iteratively. 

luese are two obvious objections to implicit schemes for our purposes: 

(i) since the interior points will be computed completely explicitly, it is 
desirable to compute the boundary points explicitly also, and (ii) computing 
normal derivatives implicitly "violates" the laws of the transmission of signals 
which are known for the hyperbolic equations under consideration. The latter is 
an objection concerned with accuracy, not stability, since purely implicit schemes 
can be unconditionally stable. 

Keeping these objections in mind, it is proper that we include implicit 
schemes in the comparison in order to understand more completely the origin and 
magnitude of inaccuracies of other schemes, as well as to be certain that no 
approach is discarded until it has had a chance to prove itself. 

2.2.5 Method of Characteristics 

Gas dynamicists and numerical analyst -e i .’.miliar with the method of 
characteristics for numerically sol vim' L -linear hyperbolic equations for 

supersonic, inviscid flow. This me ;v X on special properties of the 

governing partial differential equate. c, .w.jcn permit them to be reduced to ordin- 
ary differential equations (two independent variables) with variable coeffi- 
cients. These ordinary differential equations hold only in special directions. 
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called characteristic directions, and they reflect the facts that: (i) there a 

are preferred directions for the propagation of signals in supersonic inviscid 
flow and that (ii) these signals suffer no dispersion or dissipation (at least 
theoretically). In most instances, the method of characteristics has been the 
most accurate procedure for numerically integrating the governing differential 
equations. It is especially convenient at boundaries where other methods require 
partial derivatives normal to the boundary. It's biggest shortcoming is its 
complexity, which results from two sources. First, to employ the method only 
at boundaries while utilizing other differencing procedures at interior points 
leads to a number of interpolations that must be made in orde-. to integrate the 
equations along characteristics . Second, the equations must be combined in the 
form of characteristic compatibility equations, and this is inconvenient at best 
and time consuming at worst. The result is that boundary points require a dis- 
proportionate amount of analytical/coding effort from the analyst/programmer . 

These disadvantages are offset by the fact that the use of characteristic 

equations at the wall undoubtedly results, in general, in as accurate solutions 
. . . ** 

as xt xs possxble to achieve. That they are not more commonly used reflects 
the fact that the inconvenience to the analyst and computer coder is substantial 
and results in considerable increase in code development cost. However, if large 
numbers of production runs are to be made, these additional costs can be more 
than offset by the decreased cost per run for given accuracy that can be antici- 
pated when characteristic procedures are used at boundaries. 

Characteristics methods in one form or another have been used by many 

*** 

investigators, including Abbett, Barnwell, Coakley and Porter, Kentzer, 
Moretti, Moretti and Abbett, Moretti and Bleich, and Kutler and Lomax (1971). 

2.2.6 Miscellaneous Methods 

In this category are methods which do not naturally fall into the categor- 
ies of more standard techniques already described. Some of these procedures are 
not promising enough to receive more than passing attention; others show consid- 
erable merit. 

Extrapolation 

Probably the most simple is a simple extrapolation procedure where the 
solution at the boundary at x = x Q + Ax is obtained by extrapolating the solu- 
tion at interior points. While it is possible to construct problems in which 

— * 

From both analytical and computational points of view. 

This will be demonstrated in Section 4. 

*** 

Reported in Edelman, et al. 
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this procedure is adequate , in general it will yield completely unsatisfactory 
results, especially when the solution changes rapidly near the wall. 

Skoglund and Guy have employed am extrapolation procedure which is 
somewhat more complicated than the most simple approach. 

Kentzer* s Scheme 

A rather clever approach, first proposed by Kentzer (1970) and also used 
by Barnwell (1971), is based on the method of characteristics. The basic moti- 
vation seems to be to get most of the advantages of utilizing characteristics 
compatibility equations without some of the disadvantages . We have seen that 
there are basically two disadvantages of characteristics methods: (i) the added 

analytical work required to put the partial differential equations into charac- 
teristic form, and (ii) the additional special coding required, including a 
number of interpolations to get data at specific points on a characteristic. 

Kentzer' s approach eliminates the second disadvantage while retaining the first 

* 

and, hopefully, the accuracy of a complete characteristic procedure. 

The approach in its simplest form is best illustrated by looking at our 
basic example, supersonic, inviscid, two-dimensional flow. The governing 
differential equations for the momentum equations in characteristic form are 
are 


— dj> pu 1 2 dx 

dx ‘ ^RTT- T dx 


= 0 


( 1 ) 


along the characteristic directions 


± dy uv ± a Vu 2 v 2 - 
1 dx 


.2 - a* 


( 2 ) 


with t = v/u being the local streamline slope; Here 


d_ L + L - L + i L 

dx * cL» dx 3y “ dx A 3y 


(3) 


is a total derivative with respect to the x along the characteristic. The 
standard procedure for implementing this characteristic compatibility equation 


1 

We will see in Section 4 that it is not completely successful with regard to 
accuracy. 
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at a solid boundary is indicated in the sJ.atch. 
the characteristic between points D and 
C, where D is somewhere between the 
mesh points A and il . Si .ice the bound- 
ing geometry is specified, t ir known. 

With t d and p D obtain d by interpolat- 
ing between poii ts A and B, we have 


The equation is written along 



P C P D 


+ b(x c 


) 


Sketch 2 Method of Characteristics 
at Boundary Point 


v with b an average value of puVv'M 2 - 1 between D and C. The inconvenience and 
computational time involved in locating point D and the data there are obvious, 
particularly when there are more than two independent variables. Kentzer's 
idea is to write vhe c^mpatability equations in partial differential form at 
point A. Thus, we have 



(4) 


or 






(5) 


with 


X 


uv ± a Vu 2 + v 2 - a 2 
u 2 - a 2 


( 6 ) 


and the rest of the right-hard side evaluated at A. The interesting point is 
that, while Equation (5) hoxds generally, in the method of characteristics the 
derivatives in Equation (4) are treated as directional total derivatives in the 
directions having slopes X*. The hope is that by combining the equations in 
this special way but still retaining their partial differential character, the 
solution will be obtained with accuracies typical of method of characteristics 
calculations while retaining the simplicity of standard explicit partial deriva- 
tive calculations. 
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* 

Combined Equations 

Another approach, which has been used by Thomas and Thomas, et al., is 
to combine one or more of Euler's equations with the surface boundary condition 
Ir differential form. The object is to obtain a direct coupling of the surface 
boundary condition with the partial differential equation (s) to be integrated. 
There are many possible variations on the basic theme, and we look at some of 
the. In Section 4.3.5. 

* New Technique 

The new procedure developed in this study is analytically simple, easy to 
code, and computationally fast. It is a predictor-corrector scheme in which the 
predictor step is a regular predictor step in MacCormack's scheme. All dependent 
variables are computed in the predictor step. The corrector step is a simple 
expansion or compression wave, whichever is necessary, which turns the velocity 
vector parallel to the wall. It is discussed in more detail in Sections 4. 3. 7. 3 
and 4.6.2. 


IE 

Actually, Kentzer's scheme is one of this group, but we consider it somewhat 
separately because of its unique features. 
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SECTION 3 

WHAT TO LOOK FOR IN A COMPUTATIONAL PROCEDURE 

When confronted with the question, "What should I look for in a computa- 
tional procedure?" one's spontaneous reply is apt to be "accuracy". In this 
section we discuss the elements which determine the accuracy of a technique and 
how we can make preliminary evaluations of a technique prior to testing it. 
Moretti (1969) has pointed out the importance of considering the computational 
procedures used at boundary points as well as at interior points in evaluating 

a numerical technique. That point is certainly valid, but it does not mean that 

* 

we cannot examine and evaluate separately the components of a technique, but 
rather that in evaluating a computer code and predictions it will yield it is 
necessary to consider the contributions of all components of the technique, 
including calculation procedures used at computational boundaries as well as 
those used at interior points. Furthermore, though many of the subsequent 
arguments have much broader validity, the discussion is directed toward computa- 
tional procedures at solid wall boundary points of supersonic, inviscid flows 
where interior points are computed with MacCormack's scheme. Unless otherwise 
stated, we will be considering two dimensional, steady flow. Some discussion of 
unsteady and three dimensional flow calculations is given in Section 4.6.2. 

Obviously, a good computation procedure must be able to convey numerically 
information incorporated in the physical/mathematical problem which is being 
solved. That is why the method of characteristics is such a good procedure for 
hyperbolic problems in two independent variables. That is also why implicit 
integration procedures, though they possess nice stability characteristics, 
are not, in general, nearly as well suited to hyperbolic problems as are explicit 
procedures. The implicit procedures numerically transmit signals in directions 
quite different from those directions that the partial differential equations say 
signals ought to travel (the characteristic directions) . In supersonic inviscid 
flow, we know that the signals are transmitted along characteristics. This does 
not mean that we must use the method of characteristics to compute such a flow, 
but that any chosen method ought to maintain, in some fashion, that basic mechan- 
ism of transmitting signals. Implicit schemes do not do this. The predictor- 
corrector scheme of MacCormack evidentally does. So what does this tell us about 
numerical procedures for computing solid boundary points in supersonic, inviscid 
flow? 

l 

Which is exactly what we are doing in this study. 
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In the first place, there are two sources of information which must be 

1 

accounted for: the interior flow and its influence on the boundary, and the 

boundary geometry. We must account for the effect of each individually as well 
as their mutual interaction. That is why we can expect a rimple, completely 
explicit, integration procedure to be quite poor. Although it does have adequate 
mechanisms to account for the development of the interior flow near the wall, 
it has no mechanism to account for, during a step, the interaction between the 
wall geometry and the interior flow. Since that interaction is precisely what 
determines the wall pressure, we do not expect accurate results from such a scheme. 

Recently some investigators (c.f. Kentzer, Thomas, et al.), have attempted 
to overcome this shortcoming of explicit integration schemes by combining the wall 
boundary condition (normal velocity component = 0) into one of the partial differ- 
ential equations being integrated. The idea has merit, but, as we shall see 
in Section 4, the resulting accuracy depends very crucially on which equation 
(or equations) is chosen and how the boundary condition is invoked. 

Others (Kentzer, Barnwell) have suggested that accuracy can be improved 
by putting the partial differential equations in special forms prior to integrat- 
ing them along the wall. The idea here is that the ability of the equations to 
transmit information well depends not only on the way the differential equation 
is approximate by a difference equation but also on the form of the differential 
equation. The overall validity of this assertion is clear, but we expect to often 
find it difficult to predict a priori the relative merits of different ways of 
writing the equations. However, in general one can probably expect improved 
accuracy if the form of writing the equations somehow gives preference to direc- 
tions in which signals are transmitted (i.e., characteristic directions). 

This brings us to another interesting point, the idea of second order 
one-sided derivatives and implications with respect to domain of influence and 
domain of dependence of a point. Consider a supersonic flow going from left to 
right with a computational grid as shown in the following sketch. The solution 
is known at station i and it is to be determined at i + 1. 


Sketch 3 

Computational Grid 



Mach lines 
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We focus our attention on the problem of determining the solution at (i + l,j) 
by an explicit scheme. In a standard, second order explicit scheme (e.g.. Lax 
Wendroff, MacCormack) the solution at (i + l,j) will depend only on data at the 
three points (i,j + 1), (i,j), and (i,j - 1), which is compatable with the laws 

of propagation of signals and domains of influence and dependence. Now, suppose 
that instead of using those three points, we use the points (i,j), (i,j + 1), 
and (i,j + 2). Then, though we would be using second order accurate y derivatives, 
it is obvious that we are using information to determine (i + 1 , j ) which comes 
from outside its domain of dependence . (Note location of Mach lines in above 
sketch.) Obviously, we can think of the three point one-sided derivative as, 
in a sense, extrapolating the solution from (i,j + 1) and (i,j) to (i,j - 1) 
but we should not expect as high accuracy from the second procedure as from 
the first, and we should not be surprised if such a procedure is unstable. 

With these thoughts in mind, we can anticipate pretty well the character- 
istics of a calculation procedure even before testing it. We just need to answer 
the following question: Does the computational scheme have numerical mechanisms 

to account for the physical interactions which we know are present? If it does, 
then it should be one of the better schemes. If not, in the absence of fortui- 
tous blessings, we expect it to be in the ranks of the poorer schemes. These 
observations are not original, but they are important. Their validity is strik- 
ingly illustrated by the results given in Section 4. 
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SECTION 4 

COMPARATIVE RESULTS FOR TWO SAMPLE PROBLEMS 
4.1 GENERAL COMMENTS AND DESCRIPTION OF THE SAMPLE PROBLEMS 

Now that we have reviewed the general attributes a good computational 
procedure should have and looked at the characteristics of individual categories 
of procedures, it is time to compare the predictive abilities of the various 
procedures. The objective of the test problems is not to be exhaustive, nor is 
it to test the various methods in situations favorable to certain ones. Rather, 
it is to test each case in situations which are not necessarily favorable but 
which are representative of what is required of a technique in every day working 
codes. To do this we will consider two sample problems, a simple expansion and 
a simple compression, both for two dimensional, supersonic flow (see sketch 4). 


M = 3 

► 

Y = 1.4 
Uniform 



a) simple compression o) simple expansion 

Sketch 4 Schematic . " Sample Problem 

Expansions and compressions are, of course, the heart of supersonic flow fields, 
and these two sample problems will serve as valid tests of the relative merits 
of the various procedures. Furthermore, although these flows are comparatively 
simple, computationally they are quite challenging, as we shall soon see. 

Finally, they enable us to make the comparison against a valid standard since 
we can obtain the exact solution along the wall boundary. 

In order to present the results in as useful and clear a manner as possi- 
ble, this section is arranged as follows: First (Section 4.2) we outline the 

overall computational procedure, with some special attention paid to the boundary 
point computation. Of particular interest here are the comments relating to the 
use or non-use of the known surface entropy and total enthalpy. Next (Section 4.3) 
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we discuss the details of the surface boundary point computational procedures 

being tested. This is an important section since often small differences in 

* 

a computational procedure can sometimes significantly affect the result. Also 
in that section we introduce working abbreviations for the calculation proced- 
ures. The results for the simple compression are presented and discussed in 
Section 4.4, those for the expansion in Section 4.5. Additional comments and 
observations (Section 4.6) include a discussion of implications for three 
dimensional and time dependent flows, special comments vis e vis the envelope 
shock and its effect for the compression problem, and some discussion of the 
accuracy of the solution in the interior region. 


4.2 GENERAL OUTLINE OF THE COMPUTATIONAL PROCEDURE 

4.2.1 Basic Equations and Finite Difference Grid 

We are considering an ideal gas with constant specific heat ratio, y. The 
general flow direction is from left to right. The geometry for the two sample 
cases is shown in sketch 5. Initially, for x < 0.0, the flow is uniform and 




Sketch 5 Detailed Schematic of Sample Problems 


parallel to the x axis, with M = 3. An initial data line of 21 points stretches 
from the lower wall to the upper wall, the latter always being horizontal. 


H 

In particular, note the discussions of standard predictor -corrector procedures, 
Thomas' procedures, and Kentzer's procedure and the applicable results. 
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The basic equations for the steady two-dimensional flow of an ideal gas 
(isentropic flow, constant isentropic exponent) are 

continuity: (pu) + (pv) « 0 (1) 

x y 

1st momentum: puu + pvu + p =0 

x y x 

2nd momentum: puv + pvv + p =0 

x y y 

entropy s uS x + vS^ =0 (4) 

state: S = lnp/p Q - ylnp/p o = P - yR (5) 

For steady isentropic flow, 


( 2 ) 

(3) 


up x + vp = -=j- (up x + vo ) 


( 6 ) 


with a 2 = yp/p. Combining (6) and (1) gives the continuity equation in terms of 
derivatives c.c pressure and velocity components. 


pu x + ov y + a 7 (up x + v Py^ = 0 


(7) 


In ord^r that mesh points always fall on the lower wall, the physical 


domain is nu^ped onto a computational domain as shown in the following sketch. 

’t 

1.0 



C L / L / -X - L L-L 


0 ‘ T ~r~) — 7 7 — 7 7 7 — Tl* f 
Physical Domain Computational Domain 

Sketch 6 Physical and Computational Planes 


The computation is terminated before the effect of the upper wall can be felt. 
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The transformation is 


K ■ x 



where the lower and upper walls are described by y = B(x) and y = T(x) respec- 
tively and where 6 (x) = T (x) = B(x). Then, the equations of motion become 



Equations 8 and 9 are both linear (bu : linearly independent) combinations of 
Equations 1 and 2; thus they both contain information about the transport of x 
momentum as well as conservation of mass. Equations 10 and 11 are, respectively 
Equations 3 and 4 transformed direction to (£, n) variables where 


P = An p/p t 
R = An p/p t 
S = P - yR 

a 2 = IB 

P 


\ 

> (12) 
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For later use, we note that 



( 13 ) 


4.2.2 Interior Points 

At interior points, the equations will be integrated following the scheme 
of MacCormack as applied to equations in non-divergence form. MacCormack's is 
a predictor corrector scheme which is easily described by looking at the equation 


f C + f n • 0 


(14) 


Then, letting i be the £ index and j the n index, to obtain the solution at a 
point (i + 1 , j) using data known at a preceding £ station (i, all j) we proceed 


predictor 


f(i + 1, j) = f (i, j) ~ A£ 


[ f(i, j + 1) - f (i, j) 
An 


(15a) 


corrector 

f (i + 1, j) = i 



1 

] 

f(i, j) + f(i + 1, j) - A£ 

f (i + 1, j) - f (i + 1, j - 1) 

I 

An 

( 


- 

) 


(15b) 


For further details about the method, see the two papers by MacCormack. 

4.2.3 Lower Wall Boundary Points 

Along the surface streamline there is a certain degree of flexibility 
since the entropy and total enthalpy are known. In almost all the calculations 
we will use that information, as well as the known body slope, in order to solve 
as few partial differential equations as possible. Thus, usually we will be 
solving only one differential equacion along the wall, and that one will be an 
equation for ^ . isure. In other instances (e.g., Thomas' scheme) the basic 
scheme requires that more than one differential equation be integrated along 
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the wall. In such instances we can expect errors in surf act entropy and/or 

t 

total enthalpy to develop. We proceed in this manner because of our two-fold 
purpose: i) to compare the relative merits of various computation procedures 

currently in use, and ii) to evaluate the best- results that can be achieved 
within a particular c?.ass of schemes. For instance, take the simple explicit 
integration procedures. We could integrate all four differential equations, 
continuity, x and y momentum, and entropy. Or we could integrate one equation, 
say for pressure, and use the three relations for entropy, total enthalpy, and 
wall slope to calculate the density and two velocity components. If the differ- 
ential equation to be integrated is chosen optimally, we expect the second approach 
to yield better results since we are solving three of the four equations exactly. 

4.3 DETAILS OF THE COMPUTATION PROCEDURES USED AT SURFACE BOUNDARY POINTS 
4.3.1 General Comments 

The details of the various surface boundary condition calculation procedures 
are given in this section. In all, we will consider 25 different procedures, 
some of which are variations on various themes. 

We use simple notation. In the finite difference grid, a point is denoted 
by two indices, (i,j), where i is the index which increases in the C (abscissa) 
direction and j the index which increases in the n (ordinate) direction. At the 
initial station, i = 1. In order to invoke reflection type boundary conditions, 
the j index for lower surface mesh points is j = 2; thus j = l is reserved for 
the row of "image" points needed when reflection type calculations are being used. 

The various procedures are given abbreviations which permit a concise and 
descriptive graphical presentation of the results of Sections 4.4 and 4.5. 

In many of the schemes, we simply integrate Equation 8 to obtain the 
change in surface pressure between stations i and i + 1. In such cases, we have 

P i + 1,2 ' P i,2 + t™ 3 * < 16 > 

where RHS is some representation of the term 
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It is assumed that the reader is familiar enough with standard notation to knew 
that, for instance, in a simple explicit scheme with two-point one-sided n 
derivative?, we would have 


cP 


n 


uv n 



(18) 


etc. 

In order to calculate the surface boundary point at (i + 1,2) from known 
data at (i,j), j 2, we assume the following to be known: 

i) The complete solution at every mesh point at station i, [i.e., for 
all (i,j), j > 2] 

and, if necessary, 

ii) The complete solution of the predictor and corrector steps at all 
interior points at station i + 1 [i.e., for all (i + l,j), j >. 3J 


4.3.2 Reflection ( Image Point Procedures ) 

Reflection procedures are probably used more than any other single class 

of calculation procedures. We have considered six variations on the basic theme. 

* 

R-l R eflection in the computational plane 

All data at boundary points computed as at interior points. 
Velocity component normal to th wall not set =0. To compute 
station i + 1 from data at station i, set 


P. , = P. . 
i , 1 i , 3 


U i,l = U i,3 


v i,l = “ v i ,3 


S i,l = S i,3 


(19) 


Here R-l is the abbreviation we will use for the first reflection scheme, R-2 
for the second reflection scheme, etc. 
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R-2 Reflection in Physical Plane 

All data at boundary points computed as at interior points. 
Effect is to set, at station i, the velocity component normal to 
the wall = 0 at the wall. 


p i,l " p i,3 


S i,l “ S i,3 


L. 1 

X'l 

= u. - (cos 2 0 . - 

X f O 1 f & 

- sin e i>2 ) 

+ 2v. -sin 

0. -cos 

X r z 

e i,2 


“ v i,3 (sin2e i,2 

- COS 2 6. -) 

\ 9 

+ 2u. -sin 

X 9 3 

0. -COS 
*1“ 

e i.2 


where 0. - is the angle the tangent to the wall makes to the 
horizontal. (See sketch below.) 


y#v 



V. , = -V. . 
i»l x,3 


R-3,R-4 Reflection in Physical Plane 

Same as R-2 for pressure, but known entropy, total enthalpy, and 
surface slope used to calculate, from P i+1 the values of p i+1 2 , 
u i+l,2' and v i+l,2* Thus ' we have from Equations 12 and 13 


S s 


i+1,2 " S i,2 ^ Pi+1, 2 n i+l,2 


H i+1 , 2 “ H i,2 q i+l 


,2 " ^ 2(H i+l,2 " h i+l,2 ) 
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v i+l,2 


+ u i+l,2 


*i+l f 2 


( 21 ) 

(concluded) 


(v/u> i+l,2 - tan 9 i+l,2> 


The difference between R-3 and R-4 is that Equations (21) are imposed only after 
the corrector steps in R-3 and after both predictor and corrector steps in R-4. 


R-5 Kutler - Scheme 1 

This is one of two reflection schemes used by P. Kutler (1972). 
P f p , and u are reflected evenly, v is extrapolated. All data 
at boundary points computed as at interior points. 




R i,l ■ R i,3 


“i,l " U i,3 


v i,l " 2v i,2 ' v i,3 


v i,2 - u i,2 tan 6 i,2 


( 22 ) 


Apply v = u tan 0 after both predictor and correctc 
R-6 Kutler - Scheme 2 

Kutler' s second scheme is the same as R-4, except only the pressure 
is obtained by integrating a PDE (in this case. Equation 8). 
p, u, and v are obtained from S, H and 0. See R-3. 

R ~6A Same as R-6 except impose surface boundary condition only after 
corrector step, not after predictor step. 

R-7 Extrapolation 

All data at the image point are obtained by linear extrapolation 
from the interior. Thus 
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P i»l " 

2P i f 2 

- P. , 


R i,l = 

2R i,2 

- R. , 
1,3 


u i,l = 

2u i,2 

" U i,3 


V i,l = 

2v i ,2 

" V i/ 3 

(23) 


Pressure is obtained by integrating Equation 8 with MacCormack's 
scheme. p, u, v are obtained from the known S, H, and 0. 


4.3.3 Simple Explicit Integration Procedures - One-Sided Derivatives 

E-2 Equation 8 solved explicitly for P^ + ^ 2 {Simple Euler Scheme) 

using simple two-point, one-sided n derivatives, p, u, v obtained 
from known S, H, 0 


[e.g., P 


h = ( P i,3 “ P i,2) /An ^ 


E-3 Same as E-2 except 3 point n derivatives 


Lg., P„ = (- 3P A f2 + 4P. f3 - P ii4 )/2Anj 


4.3.4 Simple Implicit Integration Procedure - One Sided Deriviatives 


1-2 Same as E-2 except n derivatives computed implicitly as, for example. 


P n = (P i+l,3 " P i+1,2 )/An (26) 

Thus, it is necessary to iterate on the solution at (i+1,2). 
Iterate until relative change in pressure between two successive 
iterates is less than 0.0001. 


1-3 Same as 1-2 except 3 point one-sided n derivatives 

|®*g * ' P n = (" 3P i+l,2 + 4P i+l,3 ~ P i+l,4) /2An ^ (27) 

4.3.5 Standard (and not-so-Standard) MacCormack-like Predictor-Corrector Procedure s 

MacCormack's scheme is a noncentered predictor-corrector scheme in which, 
say, a forward n difference is used in the predictor and a backward n difference 
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is used in the- corrector . Now, let's determine how we should apply this approach 
at solid boundaries. 


Consider the equation 


Suppose that we know the solution at 
station i and we wish to continue it 
to i + 1. All points (i + l,j), 
j > 3 can be easily calculated by the 
standard MacCormack scheme. We are 
interested in the point (i + 1,2). 


5 <> 
4 o 
3 <> 


l 5 
o 4 
o 3 


TT7 2 7^-r7-T7T777-T7T7T77~A? ^ 

i i + 1 


Sketch 4 


The simplest approach would be to just use one-sided inward differences 
on both the predictor and corrector steps. Then, in MacCormack' s notation, we 
have a scheme we call PC-2-2. 


PC-2-2 


Predictor : 

f(i + 1,2) = f (i,2) 


Ax 


f (i , 3) - f (i , 2 ) 

Ay 


(28a) 


Corrector : 
f (i + 1,2) 



f (i,2) + f (i + 1,2) - Ax 


j~ f(i + 1 ' 3) a ~ f(i + ly2) j | 


(28b) 


where f (i + 1,2) is the regular predictor value at that interior point. For 
obvious reasons, we do not expect this result to be very accurate, irrespective 
of how the boundary condition has been incorporated into the differential equation. 

The obvious next step is to use second order accurate one-sided y deriva- 
tives. Then we have scheme PC-3-3. 


PC-3-3 


Predictor : 

Id + 1,2) = f (i,2) 


Ax 


~ -3f (i,2) + 4f (i , 3) - f (i , 4) 
2 Ay 


(29a) 


f (i + 1 


Corrector: 

,2) = | jf (i. 


2) + f(i + 1,2) - Ax 


[= 


3f(i + 1,2) + 4f(i + 1,3) - f(i + 1 

2 Ay 

(29b) 


4 


One might expect scheme PC-3-3 to be better than scheme PC-2-2. 
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To carry the investigation somewhat further, let's look at the MacCormack 
Scheme at an interior point. It is (in one form) : 

Basic Scheme 

Predictor : 

I 

Corrector : 
f (i + l,j) = | 

Comparing the basic scheme with schemes PC-2-^ and PC-3-3, we note that the cor- 
rector for scheme PC-2-2 is identical to the basic corrector. One might then 
reason that it is sensible to use the 'more accurate' predictor of scheme PC-3-3 
with the corrector of scheme PC-2-2 in trying to construct a calculation scheme 
at the wall consistent with that used at interior points. Thus, we have a 
third scheme: PC-3-2. 
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and f(i,j + 2). But what is (32)? It is 

i 

just the derivative of f with respect to 
y at the midpoint between (i,j - 1) and 
(i,j). So, to determine the term to replace 
(32), let's look at the usual one dimensional 
difference line segment (sketch 8) . The 
correspondence between points in sketches 
7 and 8 is 


4 <> 

3 o 
2 o 
1 

Sketch 8 


Sketch 7 
j + 2 

j + 1 
j 

j - 1 
j " 2 


-* — *- 


Sketch 8 
5 
4 
3 
2 
1 


Term (32) is essentially, as far as sketch 8 is concerned 



(33) 


the second order derivative of f with respect to y at the midpoint between points 
2 and 3. So, all that is necessary is to represent f^ to the same accuracy as 
(33), but ._ng points 3, 4, and 5. Following the usual procedure, we have 

*3 - *23 + *23 (¥) + I *23 (¥)‘ 

*«“*23 + *23 (!^) + I*23(! i l’) a (34 > 

*5 " *23 + *23 ( I 4 *) + I *23 (I A *) 

Eliminating ±22 and f from these three equations yields the desired equation for 
f^ in terms of f^, f^, fg. 


" 2f 3 + “ f 4 - f 5 
x 23 = Ay 


(35) 
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(Compare this with the usual three point, one-sided second order derivative used 
above in PC-3-3). Thus, we construct PCN-3-2, which we expect is the predictor 
corrector scheme for the boundary most consistent with MacCormack's scheme at 
interior points. 


PCN-3-2 


Predictor: 

I(i ♦ 1,2) = fd,2> - Ax [ ~ 2t(1 ’ 2) * ^- 3) - - f(1 - 4 >3 


Corrector : 
f (i + 1 


,2) = | jf (i. 


2) + f(i + 1,2) - Ax 


[fJL + 


1,3) ~ f (i + 1,2) 
Ay 


(36a) 


(36b) 


All four of these predictor corrector schemes have been tested. As we 
shall see in Sections 4.4 and 4.5, the results are somewhat surprising and quite 
enlightening. In all four of these "standard" predictor-corrector schemes, the 
known surface entropy, total enthalpy, and slope were used to calculate surface 
density and velocity components from the computed pressure in both predictor and 
corrector steps. Thus Equation 8 is the only partial differential equation inte- 
grated along the surface. 


4.3.6 Method of Characteristics 
MOC - Method of Characteristics 

Interpolations and solution made in physical plane. Data at the * point 
obtained by linear interpolation between (i,2) and (i , 3 ) (see sketch, 
following page) . 

X = uv_;_a Vuljf v* - a 2 (37) 


]?i+l 2 computed from campatibility equation 



in finite difference form 

P i+1,2 - P* + 5(T i+l,2 - T *> 

X and b are averaged between * point and (i + 1,2). 


(38) 


(39) 
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S i+1,2 = S i,2 * p i+l, 2 ^ h i+l f 2 

t 

H i+1,2 ' H i,2 

q i+l,2 = V 2<H i+l,2 " h i+l,2 > 1401 

(v/u) i+l,2 " ian 6 i+l,2 

v i+l,2 + “1+1,2 " '4+1,2 



4.3.7 Miscellaneous Procedures 

In this category we have those methods which do not fit into any of the 
preceding categories. Though a diverse melange, this group is the one where we 
expect the potentially greatest return as far as accuracy vs . complexity are 
concerned . 


* 

4. 3. 7.1 Kentzer's Method 

Pressure from the characteristic-like partial differential equation 



where 

3 1 3 _ 

”5y “ 7 3n 


(41) 


(42) 


The ri derivatives are two-point, one-sided derivatives explicit in x. 
That is, to go from i to i + 1, 


lE = 1 (p i/3 " p i,2 ) 

3y ST At) 


(43) 


8t _ tan 9 1+1 ~ tan 6 i 
"5x Ax 

etc. X and b are evaluated at (i,2). Then 


( 44 ) 


J 

For a predictor-corrector variation of this scheme, see Section 4.3.7. 5. 
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S i + 1,2 * S i,2 


H i+1,2 



( 45 ) 


p i+l,2> u i+l , 2 ' v i+l,2 are com P' Jt6d from 0 i+l,2 

4. 3. 7. 2 Thomas ' Procedure (see Thomas and Thomas, et al.) 

The basic idea is to replace one momentum equation with the boundary 
condition of zero velocity component normal to the lower wall. This is done 
by replacing the (x,y) velocity components, (u,v) , with components (U,V) parallel 
to and normal to the wall. Since V is everywhere normal to the wall, V = 0 
Also, 


U = u cos 8 + v sin 0 (46, 

U E = U x - | U n (47) 

Up = vi cos 0 + v sin 0 - (u sin 0 - v cos 0)0 

where 

u x " u £ - I “r, 

(48) 

v x = V C ' I v n 

Equation (47) for U is integrated alone with Equation (P) for pressure, as out- 
lined by Thomas and by Thomas, et al. 0 = d6/dx is knovrn from the specified 

velocity distrib 1 tion. V = 0, and p = complete the integration along 

the lower wall, since 


u = U cos 0 - V sin 0 


v = U sin 0 + V cos 0 


(49) 


We try two versions of this scheme, identified as Thomas-2 and Thomas-3. 

Thomas-2 - Equations 47 and 8 are integrated using two-point one-sided 
n derivatives computed explicitly ir the £ coordinate 

Thomas -3 - Same as Thomas- 3 except standard three-point one-sided ri 
derivatives are used. 
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4. 3. 7. 3 Combined Equations 


The analysis used by Thomas employs the specified variation of 0 (x) 
in computing the velocity changes by integrating the streamwise momentum equation. 
No use is made of the known variation of 6 (x) in integrating the equation for 
pressure, nor is there any direct coupling between the continuity equation and 
one or both of the momentum euqations. We know that the primary physical force 
determining surface pressure variation is momentum exchange away from the wall, 
and it is natural to consider combining the pressure equation (Equation 8) with 
one or more o* the momentum equations and with the specified variation of 0 (x) . 
This is, of course, what is done in the me J -hod of characteristics and Kentzer's 
method . 

We have investigated two such additional combinations of the differential 
equations and surface boundary conditions. They are identified by the prefix CE, 
meaning Combined Equations. 

CE-1 Inspection of Equations (3) and (7) leads to one obvious possibility. 

Multiply (7) by -v and add to (3) to get 


p(uv x - vu x ) 




= 0 


Noting that the streamline slope is 


given by 


and 


we get 


T 


V 

u 



(51) 

(52) 

(53) 


This becomes, after converting to 


(£,n) and P, 


P 






(54) 


where x r can be computed directly from the specified wall geometry. A quick 
inspection of Equation 54 gives considerable cause for worry, «tince x appears 
twice ir the denoit * nator , and x can be zero. This reflects the fact that 
Equation 54 embodies only the continuity equation and the normal momentum 
equation for determining the induced streamwise pressure gradient. The most 
important equation for determining that pressure gradient, the streamwise 
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momentum equation, was not used in obtaining Equation 54. Nevertheless, a calcu- 
lation was made for the simple expansion problem utilizing Equation 54. The 
results were quite poor, as was expected, the pressure Leing low by more than 
an order of magnitude in the first step. That large error is directly attribut- 
able to numerical inaccuracies resulting from the t in the denominator. 

Four Other CE Schemes 

Following the same basic spirit discussed above, both momentum equations 
and the continuity equation are combined, and the specified variation of 0 (x) is 
used to simplify the resulting equation and to 
changes in wall inclination and the streamwise 
multiplied by -u and added Equation 2 to get 

- - 5? Py + (l - 

Equations 55 and 53 are added to yield 

( , u 2 + uvi„ A /, v 2 + uv\_ 

1 ' — il + \ 51 — Py 


give the direct coupling between 
pressure gradient. Equation 7 is 


f) 


= 0 


(55) 


+ pU 2 (' - T ) = 0 

K x y 


(56) 


Converting to (£,n) and P we get 


P 





(u 2 + uv)c + v 2 + uv 


h 


i - 


u 


2 

~U 


+ UV 
"5 


(57) 


which is a well behaved, simple partial differential equation which can be very 
easily integrated. 

We have integrated Equatior (57) in four different ways which we identify 
as CE-2, CE-3 , CE-PC-2-2 , and CE-PC-3-2. 
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CE-2 


Equation (57) is integrated in a straight forward explicit Euler type 
scheme with two-point, one-sided n derivatives computed at station i. 
Thus, 


P(i + 1,2) = P(£ + AC, 2) = P (i ,2) + A£P^ (58) 

where P^ is given by (57) and the n derivatives are computed as, for 
example , 


p = P (i ,3) - P(i,2) 

n An 


(59) 


CE-3 The same as CE-2 except three-point one-sided n derivatives are used, 
e »g* , 


p = ~3P (1,2) + 4P(i,3) - P (i,4) 
n 2An 


(60) 


CE-PC-2-2 

This is a predictor-corrector type solution of Equation (f<7) after the 
fashion of MacCormack. First equation (57) is integrated as in CE-2 
to get the predictor value, P(i + 1,2). Then we get the n .derivatives 
in (57) by differencing (one-sided differences) the results of the predic- 
tor step at station (i + ij . Then we combine the initial value, the 
predictor value, and the corrector as in MacCormack' s scheme 

P (i + 1,2) = j jp(i,2) + P(i + 1,2) + AC-RHsJ (61) 

where RHS is the right hand side of (57) with n derivatives and coeffi- 
cients computed from the predictor solution at (i + 1) . Two-point one- 
sided n derivatives are used in both th'’; predictor and the corrector steps . 

CE-PC-3-2 

The same as CE-PC-2-2 except three-point one-sided ri derivatives used 
in predictor steps. (Two point n derivatives still used for corrector 
steps . ) 

4 . 3 . 7 . 4 Euler Predictor, Simple Wave Corrector 

E P-SWC Second order (and higher) accurate numerical procedures for inte- 
grating the equations at interior point obviously incorporate the 
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mathematical/physxcal signal propagation phenomena so well characterized 
by the theory of characteristics . The predictor/corrector procedures 
(c.f. MacCormack's Scheme) can be thought of as, in essence, computing 

the solutions for two simple waves , the solutions of which are summed to 

* 

yield the complete solution. At the wall, it is not possible to do a 

regular prei ♦•or-corrector , since there is no flow on one side of the 

** 

wall. in .ajl boundary condition replaces "the information feeding 
in from the other side". The problem is to devise an integration scheme 
incorporating these phenomena. 

Consider the simplest finite difference scheme for integrating the 
continuity and both momentum equations (p, u, v are dependent variables) . 
This scheme was utilized for the continuity equation in scheme E-2. By 
comparing 0° = v/u at x = x Q + Ax with the known surface inclination, we 
can evaluate the error in the integration. Now, we can superimpose on 
the solution a simple wave (compression or expansion) to "turn the flow" 
so that v/u will be correct. Corresponding to such a simple wave is a 
pressure increment, Ap, which is to be added to the pressure there. 

This procedure can be thought of as a predictor-corrector procedure with 
the corrector a simple wave. In the calculation, the corresponding 
pressure increment was computed from Equation 174 of NACA Report 1135 
for expansions through small angles, Av. Terms through (Av) 2 were used, 
so the equation is also valid for compression waves. The expression is 


*2 

Pi 


1 



(y + l)M!f - 4(M 2 

(Av) + yM, 

4(M 2 - l) 2 


1 ) 

(Av) 2 


(62) 


Av is positive for an expansion and negative for a compression. In the 
reported calculations, the predictor step was computed according to scheme 
E-2, except that all four of Equations (8-11) were integrated in the pre- 
dictor step (using two-point one-sided n derivatives) . Then, the pressure 
change to correct v/u was computed from Equation (62). Finally, p, 
u, and v were computed from p and the known surface S, H, and streamline 
slope. 


While this is not rigorously correct, the argument is an aid in helping to 
understand the important elements of computational procedures. 

★ * 

We omit reflection techniques from consideration here. For further comments 
on these procedures, see Section 2. 
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4. 3. 7. 5 Kentzer's Scheme: Predictor-Corrector Version 

Moretti (1971) has used a predictor-corrector version of Kentzer's 
scheme. First equations (5) and (6) are transformed to the computational plane, 
Thus 




(5) 


becomes 


3p , c + X 3p _ - / 3 t , c + X 3 t \ 
3£ 6 an “ \ 3€ $ 3n / 


(63) 


On the lower wall, c = 


-b = - — , so we have 

X u 


If - A = b (li 


+ A 


*1\ 

3n / 


(64) 


where 


^ a y(”u + 0 along the lower wall 


Now 


so that 


b = 


y/& 


a v 

pu 2 j 3 t °u 1 / \ 

and = tt- = — t ( uv - vu i 

p j 3n 3n u z V n r\) 


(65) 


35 






pu 2 3x pA 

~ ^ ./TP" 


( UV„ - vu J\ 


( 66 ) 


or 


3P 

35 


A + 

A 3n 


pu 2 3 t 


yj M 


2 - 1 


H + 71 


1L 




(UV - VU ) 

— V * V 


(67) 


In version of Kentzer's scheme, equation (67) is used in a MacCormack- 

like predictor-corrector calculation. Thus, 


and 


P(i + 1,2) = P (i,2) + P^ (i , 2) A£ 


P (i + 1 


,2) = \ |p(i. 


2) + P(i + 1,2) + P^(i + 1 


,2)A?j 


a) 

b) 


( 68 ) 
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The n derivatives at i and i + 1 are two-point one-sided derivatives between 
the initial conditions at (i,2) and (i,3) and the predictor solution at (i + 1,2) 
and (i + 1,3) respectively. In the current calculation, the term 3x/3£ was 
computed exactly, whereas in the results using Kentzer's procedure it was com- 
puted by taking a two point difference (i.e., 3r/3£ = [t (i + 1,2) - t (i,2) ]/Af; ). 

4.4 SIMPLE COMPRESSION - RESULTS AND DISCUSSION 

Our interest here is the comparative ability of the various procedures to 
predict the pressure distribution on the surface of the simple compression. The 
predicted pressures are compared wiu. each other and with the exact solution on 
Figures la-lh. The initial data (pressure and density) were specified accurate 
to four significant figures. Pertinent considerations vis a vis the envelope 
shock which forms are discussed in Section 4. 6. 1.1; it might be fruitful to 
quickly skim that discussion prior to studying the results reported in this section. 

The results for the reflection schemes R-l, R-2, R-3, and R-4 are shown 
on Figure la and for R-5, R-6, R-6A and R-7 on Figure lb. Some interesting obser- 
vations can be quickly noted. Simple reflection in the computational plane (R-l) 
is disastrous, and not surprisingly so since there is no mechanism at all for 
introducing into the flow field the fact that the surface streamline must follow 
the wall whose slope is changing. Comparing schemes R-2 and R-3, both of which 
use reflection in the physical plane, we see an example of the value of utilizing 
as much information as is available. The only difference between the two runs 
is that R-2 solves the Euler equations for all four variables, p, S, u, v, while 
R-3 solves for p only with a differential equation, obtaining S, u, and v from 
' *e known surface entropy, total enthalpy, and streamline slope. The comparison 
between R-3 and R-4 is also quite interesting. The only difference between 
these two runs is that in R-4 the known surface S, H, 8 are used in obtaining the 
predictor solution as well as the corrector solution, while in R-3 the predictor 
solution solves all four partial differential equations for p, S, u, v, the known 
surface S, H, 0 only being used in obtaining p, u, and v from p in the corrector 
step. We see that R-4 is much quicker to respond to changes in curvature, but 
that it also tends to exhibit larger overshoots. This characteristics seems to 
hold for other reflection schemes as well, at least for the compression. Kutler's 
two schemes yield essentially comparable results. There is no distinct improve- 
ment of the prediction of scheme R-6 over R-5. Again, R-6a, which utilizes 
known surface S, H, and 0 only after the corrector step, lags R-6 where the sur- 
face slope changes, but exhibits far less overshoot on the ramp section. Pro- 
cedure 1-7, which obtains all the data at the image point by linear extrapolation 

★ 

from the interior, has much larger overshoot than either R-5 or R-6, the corres- 
ponding runs with Kutler's procedures. 

I 

Recall that R-5 and R-6 (Kutler's two schemes) obtain v by linear extrapolation 

and p, p, and u by simple reflection in the computational plane. 
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We can summarize the results of the reflection schemes as follows: 

1) Reflection in the physical plane is superior to reflection in the 
computational plane. 

ii) It is preferable to obtain p by integrating a differential equation 
and p, u f v from known S, H, and 0. 

iii) Schemes which are most responsive to changes in surface curvature 
(note - not slope) (e.g., R-4, R-5, R-6, R-7) tend to exhibit 
larger overshoots than schemes '.ich are not so sensitive to changes 
in surface curvature (e.g., R- # . R-6a) . 

iv) None of the reflection schemes tested yields particularly good results. 

The results of the two simple explicit (E-2 and ’-3) and the two simple 
implicit (1-2 and 1-3) schemes also show some interesting behavior (Figure 1-C) . 

First note that the predictions of both E-3 and 1-3 using three point one-sided 
second order accurate derivatives are much worse than the corresponding predic- 
tions of schemes E-2 and 1-2 which use two-point one-sided first order accurate 
derivatives. This is an illustration of a point that is sometimes missed. Just 
because one uses higher order difference expressions it does not mean that the 
solution is more accurate. In fact, we know that three point derivatives violate 
the law of forbidden signals, and we know further that supersonic flow fields 
can have derivatives which are discontinuous in certain directions, namely 
normal to the Mach lines . Therefore . we not only should not be surprised by 
these results, we should expect them. Looking now at the two "lower order" 
schemes, E-2 and It 2, note that the implicit scheme (1-2) gives better results 
than the explicit scheme (E-2) even though 1-2 also, in a sense, violates the 
law of forbidden signals. However, it does tie in the changes in the interior 
with changes in the surface slope during one step, and its violation is "disper- 
sive" rather than "dissipative". Thus, though neither procedure is very accur- 
ate, the implicit schemes, compared to the simple explicit schemes, have a better 
mechanism for coupling the interior flow solution with the specified boundary 
condition. Finally, it should be noted that the explicit schemes are more prone 
to stability problems . Though the results of scheme E-3 dc not really look 
unstable, i is possible that the run is marginally stable or unstable. Because 
of what we know about the law of forbidden signals, it would nov oe at al 1 sur- 
prising if scheme E-3 has poorer stability limits. Thomas (AIAA 71-596) asserts 
as much with respect to his scheme (see the discussion below) using three point 
T) derivatives. 

The results for the standard predicr > ^rrector schemes are shown in 
Figure 1-d, In all four runs shown there, tne pressure was calculated by inte- 
grating Equation (8) with a MacCormack-like predictor-corrector scheme, except 
that both predictor and corrector results were computed with one-sided ri derivatives . 
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In both the predictor and corrector steps, the known surface entropy, total enthal 
enthalpy and slope were used to obtain the surface density and velocity compo- 
nents. None of the schemes is consistently superior. While scheme PC-2-2 
(two-point predictor/two-point corrector) responds to the initial expansion 
better than PC-3-3, the latter, surprisingly, has a much lower initial overshoot 
on the straight section. Both have about equivalent (but poor) ability to "home- 
in" on the correct solution along the straight section. A somewhat unexpected 
result is the poor predictive ability of scheme PCN-3-2, the scheme which is meant 
to be as consistent as possible with MacCormack's predictor-corrector scheme. 

The conclusion which one apparently must draw is that any attempt to deduce an 
"extended" (past the surface streamline ) flow field solution is not likely to 
be very successful. On the basis of these results, it is apparent that standard 
predictor-corrector schemes do not show much promise. 

The results of the method of characteristics calculation (MOC) are showr 
in Figure le, along with the results of the new predictor-corrector scheme (EP- 
SWC) which utilizes a simple explicit Euler predictor followed by a simple expan- 
sion or compression wave corrector. Both techniques yield quite accurate solutions, 
the maximum relative error in surface pressure being 0.28% for the method of 
characteristics calculation and 0.10% for the EP-SWC calculation (see discussion 
on errors below) . These small errors are achieved in calculations where the sur- 
face pressure is changing as much as 26% in one step! 

The results for both versions of Kentzer's method, Thomas-2 and CE-2 are 
shown together on Figure If since these four methods all use a partial differential 
equation formed by combining one or more of Euler's equations with the surface 
boundary condition. All four calculations employed two-point, one-sided n deriva- 
tives. Kentzer's basic procedure and Moretti's predictor-corrector version both 
predict the solution very well, the maximum relative errors being 0.56% and 0.33% 
respectively. This is indeed amazing since Kentzer's scheme integrates the 
compatability equation as a simple partial differential equation a3ong the stream - 
line (not along the characteristic) ! The results of these four calculations 
strikingly show the importance of how the differential equations and the boundary 
conditions are combined. Thomas essentially just made certain that his system of 

equations was complete and not redundant. In procedure CE-2 we have taken care to 

* 

couple the equations in a form which accounts, as well as possible, for the 
interaction between the interior flow field and the solid wall boundary. The 
superiority of CE-2 over Thomas-2 is quite clear. 

** 

The results of the two schemes Thomas 3 and CE-3 runs are shown m 
Figure lg. Obviously CE-3 is superior. In fact, Thomas-3 may be exhibiting 
some of the poor stability characteristics alluded to by Thomas (AIAA-71-596) . 


Within the constraints of the overall solution procedure and the requirement 
that we want to integrate the equation as a simple partial differential equatior 

The same as Thomas-2 and CE-2 except for three point one-sided n derivatives. 
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The results of the predictor corrector versions (CF-PC-2-2 and CE-PC-3-2) 
of the combined equations schemes (CE-2 and CE-3) are shown in Figure lh. These 
results are fairly good, scheme CE-PC-2-2 being slightly better than CE-PC-3-2. 
Both are better than CE-2 on the curved section, but not as good as CE-2 on the 
straight section. Both are far better than Thomas-2, Thomas-3, and CE-3. On 
this problem, none of the others is competitive with Kentzer's scheme. 

In order to better illustrate the comparative accuracy of the different 
procedures, on Figures 2a, 2b, and 2c we have plotted the relative error in 
computed surface pressure for representative schemes from each category. Recall 
that the initial data was specified to four significant figures, so relative 
errors less than .0002 must be considered fortuitous. Furthermore, it is impor- 
tant to look at the overall accuracy of a scheme; isolated points of significantly 
higher than average accuracy are also fortuitous. Using the method of character- 
istics as a standard, it seems reasonable to rate the schemes as follows: 


Typical Relative Error in % 


P - P 


exact 


x 100 


exact 


Rating 


< 0.5% 
0.5% to 1% 
1% to 10% 

> 10 % 


very good to excellent 
good 
fair 

unacceptable 


We see that (Fig. 2a) , for this problem, the method of characteristics, 
Kentzer's method (both versions), and the new simple Euler predictor-simple wave 
corrector scheme (EP-SWC) all rate in the very good to excellent category with 
errors mostly less than 0.3% and the average error less than 0.1%. It is interesting 
to note that all four methods appear to settle on an error trend for x/R > 0.3, 
though Kentzer's method is a little slow in homing-in on the trend. The simple 
standard predictor corrector (PC-2- 2) and reflection in the physical plane (R-4) 
yield results which are poor on the compression turn and the first part of the 
straight segment, but which improve considerably for x/R > 0.4. 

Simple explicit integration (scheme E-2) must be rated poor, while simple 
implicit integration (scheme 1-2) is in the fair to poor range (Figure 2b) . 

, * 

Of the remaining schemes in the combined equations class , Thomas ' 

scheme with two-point deriv itives (Thomas-2) is quite poor for predicting surface 

pressure, scheme CE-2 (Euler integration) yeidls fair results, and the predictor 

corrector scheme CE-PC-2-2, though slightly better on the average than CE-2, is 

also in the fair category. 


Recall that Kentzers method is also in this class. 
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Since it has yielded such good results , it is interesting to look at the 
new EP-SWC scheme calculations in a little more detail. In Figure 3 we have 
plotted the results of both the predictor and the corrector calculations. We 
see that in the first interval the predictor has no idea that the boundary 
slope is changing, so the entire change in the solution along the surface stream- 
line is taken up by the simple wave corrector step. Subsequently, the non- 
uniform flow above the wall carries with it information about the changing wall 
slope, and the predictor contributes a large part of the change until the straight 
segment is reached. For the first interval thereafter the predictor acts as though 
the surface slope were going to continue to increase, and the corrector brings 
in again (as in each step) the interaction between the surface and the interior 
flow. Finally, for 0.242 <. x/R < 0.34, the predictor computes an expansion which 
must be corrected with a simple wave compression. 

It is also interesting to note that in Moretti's predictor-corrector version 
of Kentzer's scheme, except for one interval, the contribution of the corrector is 
quite small. On the curved section the change from predictor to corrector is 
between 0.13% and 0.85% while the solution is changing between 14% and 26% in 
one interval. It is only in the interval where the surface curvature is dis- 
continuous that the corrector contributes a major portion of the change. In 
that interval the charge between predictor and corrector is 8.3%. Thus, the 
strength of Kentzer's approach appears to primarily reside in the form of the 
differential equation, rather than refinements in the differencing procedure. 

4.5 SIMPLE EXPANSION - RESULTS AND DISCUSSION 

The simple expansion is the second of this pair of elemental problems 
used to test these computational procedures. The initial rnesh spacing and flow 
conditions are the same as for the compression. 

The predictions of the various reflection schemes are compared with the 
exact solution in Figure 2a and 2b. Again, simple reflection in the computa- 
tional plane is disastrous. These results differ somewhat from those for the 
compression in the following respects: 

i) Scheme R-2, reflection in the physical plane with p, S, u, v all 
computed by integrating Euler's equations, shows less ability to 
home-in on the solution on the straight segment. 

ii) Scheme R-4 is again superior to R-3 on the curved surface, but it 
is now equal to R-3 on the straight segment. 

iii) The percent overshoot of R-3 is higher (19.1% vs 8.4% for the 
compression) . 

■?v) Compared to R-3 and R-4, Ku bier's schemes R-5 and R-6 do better 
than they did on the compression. 

v) R-7 again shows comparatively large overshoot. 
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The results of the two simple explicit (E-2 and E-3) and two simple implicit 
(1-2 and 1-3) schemes (Figure 4c) are a bit different than they were for the 
simple compression. On this problem there is not clear superiority of two point 
n derivatives compared to three point derivatives. In fact, overall probably 
1-3 is a little superior to 1-2 and E-3 to E-2. Both E-3 and 1-3 respond quicker 
than the corresponding E-2 and 1-2 to the initial curvature change and to the 
change vrhere the straight segment begins, but again E-3 and 1-3 exhibit greater 
overshoot on the straight segment. As for the compression, none cf these schemes 
yield good results. 

The four standard predictor-corrector schemes yield comparable results 
(Figure 4d) . Again, overall PC-2-2 seems to be somewhat better than the others, 
primarily bev.ause of its superior ability to home-in on the solution on the 
straight segment. Again, none of these results are acceptable, with errors as 
high as 26% for scheme PC-3-3. 

The method of characteristics (MOC) and the new simple-Euler-predictor/ 
simple-wave-corrector (EP-SWC) schemes yield excellent results (Figure 4c) . The 
maximum relative errors are 0.58% for the method of characteristics and 0.43% for 
the EP-SWC scheme. This high accuracy is achieved when the exact solution changes 
as much as 22% in one stepl 

The results of Thomas procedure and of the basic combined equations approach 
with Euler integration (CE-2) , both using two-point one-sided r\ derivatives, are 
compared with Kentzer's predictions and the exact solution on Figure 4f. Again 
we see that scheme CE-2 yields far better predictions for surface pressure than 
does the Thomas-2 scheme, and Kentzer's methods are still the best. Kentzer's basic 
procedure, while still pretty good, does not do as well on the expansion (Fig. 4 f ) 
pansion (Figure 4 f ) as it did on the compression (Figure If) . While doing pretty 
well on the curved section, it has about a 10% overshoot on the straight segment, 
after which it very quickly homes -in on the correct solution. Even on the curved 
section, typical errors are considerably higher for the expansion (1-7%, Figure 
5a) than they were for the compression (0. 8-1.0%, Figure 2a). There is no ob- 
vious explanation for the comparatively poorer ability to predict expansions 
than compressions. On this problem, Moretti's predictor-corrector version model 
is considerably better than the basic Kentzer scheme. The overshoot is down to 
3.4%, and on the curved section the errors are 1/2 to 1 order of magnitude smaller. 
Again, except in the one interval where the surface curvature is discontinuous, 
the corrector changes the predictor by less than 1%. Thus, we can attribute some 
of the differences between Kentzer's basic scheme, as used here, and Moretti's 
predictor to the different treatment of the derivative of the surface slope (see 
Section 4. 3. 7. 5). Additional studies of these two approaches and the treatment 
of the surface slope derivative will be reported in the future. 
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Schemes CE-3 and Thomas-3, both using three-point one-sided n derivatives, 
are compared with the exact solution in Figure 4g. The results here for Thomas-3 
are considerably better than they were for the compression, but it again shows 
a strong tendency to overshoot on the straight segment. 

Finally, the results of the two corresponding predictor-corrector schemes, 

CE-PC-2-2 and CE-PC-3-2 are shown on Figure 4h. Scheme CE-PC-3-2 is slightly 

better than CE-PC-2-2, and on this problem both are better than any of the oth^r 

* 

schemes in this category, except Kentzer's scheme. 

The relative errors of some of the schemes are compared with ec ther 
in Figures 5a, 5b, and 5c. The method of characteristics (MOC) and th jw Euler 
predictor, simple-wave corrector (EP-SWC) scheme are far superior to any of the 
others, with Kerttzer's scheme considerably better than those remaining. 

As with the expansion wa\e, the path of the EP-SWC predictor corrector 
scheme is given in Figure 6 (see Figure 3 and its discussion) . 

4.6 FURTHER COMMENTS AND OBSERVATIONS 

4.6.1 Accuracy of the Solution at Interior Points 

4. 6. 1.1 Simple Compression - The Envelope Shock 

For the compression problem it is important to consider the influence 
that the envelope shock which will form may have on the results. No provision 
was taken to treat that shock as a sharp discontinuity, so it is possible that, 
particularly with those schemes using three-point one-sided differences, the 
solution at the lower wall is affected by numerical inaccuracies resulting from 
the shock smearing. 

As a first step to evaluate this effect, a full method of characteristics 
calculation was generated with a code which computes shock waves as sharp discon- 
tinuities. The code (referred to as SUPER) , detects envelope shocks by the 
crossing of characteristics of the same family. It is a second order accurate 
solution since both characteristic slopes and coefficients in the equations are 
averaged in a solution which iterates each step. In orde~ to obtain mesh point 
spacing on the lower wall at about the same interval as in the "finite difference" 
calculation, the mesh points on the initial data line were s^ : jed at distance 
Ay/R = 0.01 apart. The calculation proceeds from mesh points on the initial 
data line (x = C.O) along right-running characteristics, and it was continued 
until the code halted because the number of points on the characteristic line 
exceeded 50, the maximum allowed. The last computed right-running characteristic 
originated at x/R = 0.0, y/R = 1.28. The surface pressure computed with the 
full method of characteristics calculation agrees very well with the exact 

J . 

Gas dynamic equations combined with surface boundary condition in differential 
form. 
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tropic compression solution (Figure 7) . The resulting flow field and some inter- 
esting additional information are shown in Figure 8. The shock was first detected 
at x/R = 0.254, at which poinr a sharp shock solution yielded a pressure ratio of 
p = P 2 /Pi = 1.129. In order to determine as closely as possible the actual abscis- 
sa of the beginning of the shock, the pressure ratio across the shock was plotted 
as a function of x/R and the plot was extrapolated back to the abscissa where 
p = P 2 /Pi = 1*0 (Figure 9)# yie^ing the initiation of the shock at x/R = 0.24. 

The right running characteristic passing through that point hits the lower wall 
at x/R =0.38. On the last computed right running characteristic, the shock was 
located at x/R = 0.379, y/R = 1.146, and the shock slope has almost reached its 
asymptotic value (see Figure 8) . Also shown on Figure 8 are the initial finite 

difference arid for the calculations reported herein, the finite difference grid 

* 

at x/R = 0.255, the left running idach Line which originates at the beginning of 
the compression turn, and the asymptotic shock angle. Note that, including the 
lower wall point, there are four mesh points between the lower wall and the shock 
in the region x/R ^ 0. 7 ~ Further downstream, there will be more points between 
the shoe* and the low -.11 (see Figure 8). 

In a study of one-dimensional time-dependent flows with shock waves com- 
puted by shock smearing, Gary compared shock speeds predicted by Lax-Wendroff 
differencing of tl Euler equations in divergence form with predictions using 
the equations in non-divergence form (as done in this study) . C *v obtained 

much g-eater errors in computed shock speed utilizing non-divergence form than 

* * 

when the divergence form was used. Kutler anticipated similar behavior in the 
present case, expecting the smeared shock location to differ signif icrntly from 
the correct value. To investigate this behavior as well as to study the effect 
of the shock smearing on the present results, the computed pressure profiles at 
three sections are plotted for Case EP-SWC in Figure 10. In Figure 10a we see 
that when the shock is just beginning to ferm (x/R - .25) the region of strong 
pressure gradients as predicted by the finite difference grid somewhat lags 
that predicted by the full method of characteristics calculation. The finite 
difference solution overshoots by about 9% the pressure profile on the compression 
side. At x/R = C . 3664 'Figure 10b) the shock is well established, and the smeared 
shock is displaced even more from the exact solution, and the pressure overshoot 
is now about 27%. .btcause of the mesh poi..t limitation, we have no complete method 
method of cbaractej. x . tics solution for all y beyond x/R = 0.37. However, since 
the computed shock angle is Try close to the asymptotic value, we can extrapolate 
it to obtain the shock location at another abscissa, for example f s: x/R = 0.52 '>4 
(Figure 10c) . Then it is aae to * onstruot the exact pressure profile and to 

"Tt 

Recall that, at a given x, th re are twenty equal intervals between the lower 

wall and the ordinate y/R = 1.4. 

*•* 

Private communication. 


42 



compare it with the prediction using the finite difference formulation with shock 
smearing (Figurfe 10c) . Again, we see the displaced shock location and pressure 
overshoot. (still 27%) predicted with shock smearing. 

Now, what is the importance of these results vis a vis this study? One 
might remark, for instance, that if the envelope shock were treated as a sharp 
discontinuity, the problem of inaccuracies due to displaced shocks and pressure 
overshoots would disappear. However, that is misleading, for the inaccuracies in 
boundary point calculations do not primarily originate with the poor shock des- 
cription. Also, no matter how the shock is treated, we must always live with 
some errors which will result in the interior calculation even if the boundary 
solutions are exact. Furthermore, we must fact the situation depicted in Figure 
10a, where the shock, as a shock, contributes a negligible part of the compres- 
sion, but where there still is the usual pressure overshoot which occurs even 
when the equations are cast in divergence form. Thus, it is just as erroneous 
to replace that smooth compression with a shock as it is to represent the compres- 
sion through a shock with shock smearing. Even if we attempt to devise a code 
in which shocks are essentially always treated as sharp discontinuities, we must 
realize that if the flow is at all complicated, weak secondary compressions and 
shocks may exist which are not accounted for within the code logic. Our task is 
to insure that in such situations relatively minor inaccuracies in the interior 
do not lead to large inaccuracies at boundaries. 

Finally the purpose of this study is to evaluate various computation 
procedures in situations they are expected to handle as a matter of course. All 
the computations are identical except for the boundary point computations, and 
the tests have been of situations typical of those faced by everyday working codes. 

It is clear that schemes using three-point one-sided derivatives will be 
much more susceptible to the small or large oscillations shown in Figures lOa-lOc. 
However, the disadvantage of three-point differences goes deeper than that. Even 
if the interior solutions were exact, these procedures will generally exhibit 
the comparatively almost undamped overshoot discussed in Section 4.4 and 4.5. 

4. 6. 1.2 Simple Expansion 

To further clarify some aspects of the interaction between solutions at 
boundary poin J 3 with those at interior points, we consider ti e simple expansion 
problem. Since the solution is just a simple Prandtl-Meyer turn, we can calcu- 
late the exact solution at any point in the flow field. First we consider one 
of the better schemes for computing boundary points. The solution using the 
EP-SWC procedure at boundary points is compared with the full method of character- 
istics solution and with the exact solution, all at x/R = 0.2648, in Figure 11. 

The agreement is very good between all three, the full method of characteristics 
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solution showing excellent agreement with the exact simple wave solution. The 
maximum error at an interior point for the calculation using EP-SWC is about 2 
to 2*1/2%, which is somewhat higher than the errors on the boundary point. 

Further downstream, the pressure profile has a sharp discontinuity in 
slope. The results of the EP-SWC calculation are compared with the exact solution 
at x/R = 0.5095, 0.6654 in Figures 12a and 126 respectively. Note that the agree- 
ment in the expansion Atjion is quite good, but that there is a tendency to de- 
velop small wiggles near the discontinuity in profile slope. However, as was 
demonstrated in Sections 4.4 and 4.5, there are schemes for computing boundary 
points which are insensitive to these small wiggles, and the wiggles need not 
have a deleterious affect on the boundary calculations. 

It is also interesting to compare the pressure profile for a couple of 
the poorer schemes with the exact solution. The results for schemes CE-PC-2-2 
and Thomas-2 at x/R = .43 are show-i - Figure 13. Note that the large errors in 
computed surface pressure do not ree ly have much effect on the interior flow, ac 
least in this case. However, it is important to remember that surface pressure 
is generally the single most important result of these calculations! 

4.6.2 Three Dimensional and Time Dependent Flows 

We can expect most of the results and conclusions to carry over to time 
dependent as well as tin dimensional steady flows. The details of some of the 
procedures will vary, of course, but the general principles we have been discussing 
and testing are not restricted to two independent variables. The advantage of the 
new predictor-corrector scheme for three or more independent variables is obvious 
since the implementation of either Kentzer's scheme or the method of characteristics 
is quite cumbersome. There is a certain degree of arbitrariness in applying the 
scheme tc three or more independent variables, but comparisons with method of 
characteristics calculations for shuttle type vehicles have been very favorable 
(Rakich & Kutler, 1972) . 

The extension to three dimensional steady flow can be easily sketched. 
Consider, for example, a cartesian coordinate system (x, y, z) with velocity compo- 
nents (u,^,w) and unit vectors (i,],£). Let the flow be supersonic in the general 
direction of the x axis (i.e., *i must be supersonic, ** and w will generally be sub- 
sonic, though they also could be supersonic) . If the surface geometry ii> given by 
an equation of the form 

F(x,y,z) = 0 (63) 
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then the unit normal to the surface is given by 


- V + V + F z* 

n = *— = n.i + n-*i + n n k (64) 

Finally, we denote the velocity vector by q. 


q = ui + vj + wk (65) 

As in the two dimensional case, we first generate a predictor step solution 
for p , p , u , rT , w using MacCormack's scheme with inward differences. Then, 

O 0.0 0 

since the boundary condition is that 


> 

q 



( 66 ) 


along the surface, we find the angle between the predicted velocity vector 


q=ui+vj+wk and the unit normal. 
^ o o o 


Av = sin” 1 



n 


(67) 


Then we impress a simple wave to turn the velocity vector through the angle Av. 
Again, associated with that wave is a pressure change given by Equation (62) . 
Next, we find the corrected velocity components from the expression 


<*1 = % 



( 68 ) 


If v.e know exactly the surface total enthalpy and entropy at the point, we can 
find the exact value of the density and velocity modulus corresponding to the 
corrected pressure, p^. In such a case, we can make a final small correction to 
the velocity vector by maintaining the direction given in (68) and scaling its 
modulus to the known value. 

It should be clear that this procedure is not completely rigorous. For 
one thing, there is no simple wave in three dimensional flow, so the use of 
equations (67), (62) and (68) is certainly an approximation. How good an 
approximation it is can only be determined after the fact. Preliminary results for 
sharp cones and Hunted shuttle vehicles at incidence have been very good. This 
should not be too surprising since one can often treat three dimensional flows 
as locally two-dimensional with crossflow effects accounted for in an "inhomogene- 
ous" term. 
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Obviously we could construct the corresponding scheme for unsteady flow 
computations. In fact, in a sense the corrector is similar to Godunov's method 
for time dependent flows. 
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SECTION 5 

CONCLUDING REMARKS 

It need hardly be mentioned that the range of computational procedures 
considered in this study is not completely comprehensive. One can, of course, 
think of quite a number of variations on schemes investigated here, and it is 
very likely that one could invent a scheme quite different to any of those we 
have examined. However, this study is complete, as far as procedures are concerned, 
in that all currently used procedures are probably represented in one of the cate- 
gories, and the schemes tested are representative of the best results that can 
be achieved in each category. 

The procedures have been compared for an ideal gas at Mach = 3 in two- 
dimensional flow. The results would be essentially the same for an equilibrium 
gas. Although we anticipate quantitative differences at other Mach numbers, 
qualitatively these results should carry over as Mach number varies. Similar 
comments should apply if the geometry of the test problems were altered. 

There are a number of observations to be made abcut and conclusions to 
be drawn from the results presented in Sections 4 . 4 and 4.5. It should be evident 
that any purely explicit integration procedure not incorporating explicitly the 
wall geometry changes possesses an inherent lag in that it has no mechanism for 
introducing into the solution for one step changes in wall geometry occuring during 
that step. Thus, we can expect such procedures to give poor results any time well 
curvature effects are important. On the other hand, we have seen that purely ex- 
plicit reflection procedures tend to recover fairly quickly once the surface geometry 
becomes straight. 

Higher order one-sided differences do not generally yield better results 
(c.f., -esults for schemes E-3, 1-3, PC-3-3, Thomas-3, CE-3) . They tend to be more 
quickly responsive to changes in surface slope, but they also show strong tenden- 
cies to large, pool / damped oscillations. As already noted, this should not be 
surprii?’ the use of three point differences violates the law of forbidden 

signa. 

' Nation procedures have given markedly better results 

th« • •£ '• .’e best procedures a ,- e the method of characteristics, 

the - • - ; . ' -a. cor simple wave corrector procedure, and the two ver- 

si ri o *■ ' i Though very competitive on the compression, Kentzer's 

meth . .. . •. ■ comparatively on the expansion problem. 
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The success of the new predictor-corrector procedure (EP-SWC) is probably 
to be explained as follows: The simple Euler predictor step accounts for the 

development of the interior flow between the surface and the next point in. The 
simple wave corrector step then accounts for the interaction of the predictor 
solution with the surface boundary. We have seen that the scheme, though quite 
simple, yields very accurate results. It has been extended to three-dimensional 
steady flow* and the results agree very well with those achieved using the 
method of characteristics in three dimensions. It requires very little comput- 
ing time and, because it is so simple, is unlikely to result in numerical prob- 
lems in difficult situations. While no formal analysis of its accuracy has 
been made , the results achieved so far indicate that it is accurate to second 
order in mesh spacing. 

Other than the developme t of the new scheme, the most important result of 
this study is the very striking demonstration of the fact that accuracy depends 
very critically on the manner in which the equations are solved. This fact is 
clearly demonstrated by comparing the results of methods E-2, Thomas-2, CE-2 and 
Kentzer, all of which used two-point one-sided n derivatives. The latter two 
procedures give considerably better results than the former two. The reason for 
the differences lies solely in how the equations are written, including whether 
and how the surface boundary condition is incorporated in the partial differential 
equation for pressure. 


The scheme outlined in Section 4.5.2 has been incorporated in Kutler s co e 
(Kutler, et al., 1972) and comparisons have been made vith a full three- 
dimensional method of characteristics calculation for a shuttle type vehicle 
at M = 7 , and 5° angle of attack (Rakich and Kutler, 1972). 
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